Models of molecular evolution

Models of molecular evolution have several components that each consider different aspects of the evolutionary process. The models that have been used to model protein evolution in studies of the deep phylogeny have largely differed according to three components: the exchangeability matrix \(E\), the rate of evolution, and the amino acid equilibrium frequencies \(\Pi\). A detailed description of each component is in Supplementary Results. The exchangeability matrix \(E\) describes the rate at which one amino acid changes to another. Exchangeability matrices used in the studies under consideration here include Poisson (all rates equal; also known as F81), WAG and LG (fixed rate matrices derived from large databases of protein sequences), and GTR (rates inferred from the data at hand). The rate of exchange can be further scaled, for example with a Gamma parameter that accommodate rate variation across sites. The approach to amino acid equilibrium frequencies includes using the frequencies observed in the data, estimating a single set of frequencies that are optimized over all sites for all taxa, or using heterogeneous models where frequencies are allowed to vary. Heterogeneous equilibrium approaches include CAT³, which is implemented in the software PhyloBayes and has been widely applied to questions of deep animal phylogeny, and the CAT-like C10 – C60 models (C models) employed in maximum likelihood (ML) inference until now has not been applied to the problem.

Models of molecular evolution are assembled by selecting different options for all these different components. The models that are applied in practice are heavily influenced by computational costs, as well as convention. For example, on the questions considered here, CAT site-heterogeneous models of equilibrium frequency have only been used in combination with Poisson and GTR exchangeability matrices. LG and WAG exchangeability matrices have only been used with site homogeneous estimates of equilibrium frequency. This is further confused by the abbreviations that are used for models. Papers often discuss CAT and WAG models as if they are exclusive, but these particular terms apply to non-exclusive model components– CAT refers to equilibrium frequency variation across sites and WAG a particular exchangeability matrix. In this literature, CAT is generally shorthand for Poisson+CAT and WAG is shorthand for WAG+homogeneous equilibrium frequency estimation. One could, though, run a WAG+CAT model, although the computation time is much longer than Poisson+CAT model. To avoid confusion on this point, we always specify the exchangeability matrix first, followed by modifiers that describe the accommodation of heterogeneity in equilibrium frequencies (e.g., CAT). Moreover, the Gamma-rate heterogeinity (G4) is used in almost every analysis conducted here. If there are no modifiers, then it is implied that site homogeneous models are used.

Gene sampling

High-throughput sequencing allows investigators interested in the phylogeny of early diverging animals to assemble matrices with hundreds or thousands of genes (Table 1). Surprisingly, gene overlap between studies has rarely been assessed. Though there are many genes with homologs present in most known animal genomes (e.g., animal BUSCO genes⁴), studies of animal phylogeny often relied on a wide variety of different bioinformatic approaches to identifying and selecting genes for their matrices. These take different approaches to sequence comparison, filtering, and orthology assessment, and gene selection can be impacted by which datasets are considered (see supplementary text for more details). As a result, the genes selected for analysis can vary widely between studies.

Outgroup and ingroup taxon sampling

Outgroup selection (the non-animal taxa used to root the animal tree) can affect the phylogenetic inference of animal evolution (Fig. 1C). Several studies have progressively removed more distantly related outgroup taxa to test the effect of outgroup selection to ingroup topology^5,6. Three progressively more inclusive clades have often been investigated: Choanomalia (animals plus most closely related Choanoflagellata), Holoza (Choanimalia plus more distantly related holozoans) and Opisthokonta (Holozoa + Fungi). It has been suggested the inclusion of more distant related (e.g. Fungi) outgroup taxa tend to recover Ctenophora-sister, whereas increasing support of Porifera-sister often recovered when more closely related outgroup taxa were used with site-heterogeneous Poisson-CAT models⁶. In contrast of outgroup selection, sensitivity to ingroup sampling has received less attention.

Data-recoding

Compositional heterogeneity (clade-specific and site-specific differences in amino acid frequencies) have been a concern in the study of phylogenetics, and has been of particular concern for the deep animal phylogeny⁷. Variation in evolutionary rates across taxa can lead to compositional heterogeneity that may lead to long branch attraction of sponges or ctenophores⁸. Feuda et al. 2017⁹ has proposed to use data-recoding method to reduce compositional heterogeneity in matrices. By discarding information on which specific amino acid is found at each site in each species and instead focusing on which groups those amino acids belong to, they sought to reduce the impact of differences between species in which particular amino acids within those groups are most frequent. By reanalyzing Chang et al. 2015¹⁰ and Whelan et al. 2015¹¹, they report that posterior predictive analyses indicate recoded analyses have better model adequacy than non-recoding analyses, and that Porifera-sister was strongly favored under all recoding strategies using GTR-CAT models⁹. However, how the recoding methods impact phylogenetic inference has not been well explored. Moreover, a recent simulation study has suggested that non-recoding approaches significantly outperformed recoding strategies¹². Thus, it is important to reevalute the performance of the impact of data-recoding compared to non-recoding methods in animal phylogeny.

Matrix taxon composition

Taxa sampling plays a critical role in phylogenetic inference, especially for the taxa whose phylogenetic position is poorly understood. However, several data matrices have relatively skewed taxon sampling, with more than 50% taxa sampled in Bilateria compared to other major animal lineages (Extended Data Fig. 1). These skewed data matrices include two studies constructed from ESTs (Dunn2008² and Hejnol2009¹³) and one study from whole-genome data (Borrowiec2015¹⁴). Importantly, several earlier studies before 2015 only include less than five ctenophore species (Extended Data Fig. 1). The rest data matrices have relatively balanced taxon composition, although the total number of sampled taxon varies. In general, increasing taxon sampling is broadly accepted to improve phylogenetic inference and aid in the placement of the earliest diverging animal lineage. Thus, it should be noted that Whelan2017¹⁵ (has the most Ctenophore taxa) and Simion2017¹⁶ (contains the most taxon sampling) have much broader taxon sampling than other matrices.

Fig. 2. Matrix taxa and gene sampling. Each of the primary matrices considered here, color coded by taxon sampling. Horizontal size is proportional to the number of genes sampled, vertical size to the number of taxa sampled.

Matrix gene composition and overlap

We first used the Metazoa BUSCO genes⁴ to evaluate gene composition from all data matrices and the result is unexpected (Extended Data Fig. 2A). The animal BUSCO dataset contains 978 single-copy orthologs (identified from 64 animal genomes) that are most likely shared by most animal genomes, and are commonly used to assess the completeness of genome or transcriptomic assemblies⁴. The BUSCO genes also are frequently used as markers in many large-scale phylogenomic studies (e.g., yeasts¹⁷, spiders¹⁸). Thus, it is reasonable to expect that these conserved genes should present at a higher percentage than non-BUSCO genes in most matrices. However, we found that gene partitions from all matrices comprised less than 50% of BUSCO genes, even with several data matrices contain less than 20 percent (Extended Data Fig. 2B). Moreover, we identified number of ribosomal RNA genes in each data matrix (Extended Data Fig. 2C) and found that Philippe2009, Nosenko2013_ribo_14615 and Chang2015 matrices possessed much higher content of ribosomal genes than other matrices (Extended Data Fig. 2D). We also annotated each gene partition against SwissProt¹⁹ and Gene Ontology databases, and results are shown in Supplementary Table 1.

Previous analyses have mainly focused on outgroup selection and/or model fits in different data matrices, and the matrix overlap across different studies has been rarely investigated since all these studies rely on their own method for orthology inference. Our initial intent of comparing matrix overlap across matrices was to construct a new matrix with shared taxon and gene sampling to test specific hypotheses about differences in support. However, our results show that surprisingly little overlap of gene samplings across all the data matrices from different studies (Fig. 2). Thus, these results strongly suggested that gene sampling largely varies across studies, indicating that different regions of the genomes are sampled with little overlap among different data matrices. Moreover, the little overlap among studies may potentially lead to incongruent phylogenetic results and difficult to directly compare results across studies.

Support for Porifera-sister and Ctenophora-sister

Fig. 3A. A summary of phylogenomic results from the previous studies and reanalyses conducted in this study. (A) A total of 135 analyses were transcribed from the literature (Table XX). (B) A total of 98 new phylogenomic analyses were conducted in this study (Supplementary Table XX).

Here we summarized all the 163 previous phylogenetic analyses from 15 studies (Fig. 3A and Supplementary Table 2). Among these studies, the main conclusions of five studies strongly in favor of Porifera-sister and 10 in support of Ctenophora-sister (Table 1). Moreover, three studies are secondary studies without any new data matrices. Interestingly, we found that all analyses strongly support of Ctenophora-sister once site-homogeneous models were used, no matter what outgroup choice, phylogenetic inference and whether data-recoding are used. The support for Porifera-sister increases with the exclusion of more distantly related outgroups and with the use of site-heterogeneious CAT (Poisson+CAT or GTR+CAT) models in several matrices, including Philippe2009, Ryan2013_est and Whelan2015. It should be noted that many BI analyses with CAT model are not well resolved (either not reach convergence or do not provide significant support for the key node of animal evolution) mainly due to computational limits of CAT models employed in PhyloBayes. The strong support of Porifera-sister can only be recovered from analyses using Poisson+CAT and recoding methods with GTR+CAT models in several matrices. Although GTR+CAT model generally has better model adequacy than Poisson+CAT or site-homogeneous models from previous studies^6,9, Porifera-sister is not always supported since several analyses strongly in favor Ctenophora-sister even with only choanofalgellates are used as outgroups (Fig. 3A). Thus, our summary of previous phylogenetic analyses strongly indicated that the trend of the support of Porifera-sister is generally increasing with the use of more complex CAT models and more closely related outgroups in several representative matrices, whereas Ctenophora-sister is supported in most other conditions.

Comparison of IQtree and PhyloBayes results

Site heterogeneity in equilibrium frequency has been a major concern in tests of Ctenophora-sister and Porifera-sister. This has been addressed with CAT models. IQtree provides a new family of C60 models that also address site heterogeneity. Given the extensive computational cost and concerns about overparameterization of CAT models²⁴, we compared IQtree C60 results to CAT results for a subset of matrices to see if they give consistent results. This would be of technical interest because it would reduce the cost of accommodating compositional heterogeneity in future analyses. It should be noted that most Phylobayes runs were converged, although several large matrices have not reached convergence after at least a month’s computational time.

Interestingly, we found strong supports for both hypotheses when using Poisson-CAT model in analyzed Choanimalia matrices. Similar to ML analyses, we found no strong support for either hypothesis in Moroz2013_3d and Nosenko2013 matrices. Consistent with previous study^25,6, we recovered strong support of Porifera-sister in matrices of Philippe2009, Ryan2013_est and Whelan2015. Interestingly, we also found a strong support of Porifera-sister in two Whelan2017 matrices¹⁵ and a weaker support of subsampled matrix of Simion2017¹⁶. The Ctenophora-sister were strongly supported in all other matrices. For the matrices have conflicted results between Poisson-CAT and models employed in IQtree, we further ran Poisson+CAT models with inclusion of more distantly related outgroups to evaluate how outgroup choice may affect phylogenetic inference under CAT model. We found a weak support of Porifera-sister in all Holozoa matrices and strong support of Ctenophora-sister in all analyzed Opisthokonta matrices. These results suggest that Porifera-sister can only be recovered in certain gene sampling in several representative matrices with the exclusive of more distantly related outgroups when analyzed using Poisson-CAT. Thus, our results from reanalyses are largely similar to previous studies and indicate that both model choice and outgroup sampling play a major role in incongruent results of animal origin in several representative data matrices.

Phylogenetic signal

The distribution of phylogenetic signal has never been quantified in site-heterogeneous models and compared with site-homogeneous models. To further explore the effect of model selection and outgroup choices in animal origin, we next quantify the support of phylogenetic signal over two alternative hypotheses (T1: Ctenophora-sister (Fig. 1A); T2: Porifera-sister (Fig. 1B)) to three representative data matrices with different outgroup choice and models in both ML and BI framework (Fig. 5). By calculating gene-wise log-likelihood scores between T1 and T2 for every gene (\(\Delta\)GLS) or site (\(\Delta\)SLS) in each matrix, we found that Ctenophora-sister had the higher proportions of supporting genes in every analysis when using site-homogeneous models (LG or WAG in both IQtree or PhyloBayes). Moreover, the outgroup choice has little impact on the distribution of the support of phylogenetic signal in sifte-homogeneous models. This finding is largely consistent with the results from past study that majority of phylogenetic signal is strongly favored in Ctenophora-sister hypothesis with site-homogeneous models in other data matrices²⁶.

In contrast, the phylogenetic signal of site-heterogenous models changed dramatically compared to the site-homogenous models. Although Ctenophora-sister also had the higher proportions of supporting genes with C60 models, we found that the phylogenetic signal decreases in many genes using C60 models compared to site-homogeneous models, especially in matrices from Ryan2013_est that nearly 30% of genes changed from strong to week signal ($lnL<2) (Fig. 5B; Supplementary Table 4). It is not clear why the phylogenetic signal decreases in many genes with a better fitting model as suggested by model testing in IQtree, but our results indicated the C60 models lack of phylogenetic signal compared to site-homogeneous models in the node of the root position of animal phylogeny. In contrast to the C60 models, the phylogenetic signal largely increased in Poisson-CAT models and more importantly, the outgroup choice has a major effect of the distribution of phylogenetic signal (Figs. 5A, C). We found that support for Ctenophora-sister largely decreases and Porifera-sister increases (~ 20% of sampled gene partitions) when exclusive of more distantly related outgroups with the Poisson-CAT models (Supplementary Table 4). The results here are further corroborated by our phylogenomic analyses in BI, which the support of Porifera-sister generally increases in Holozoa and Choanimalia matrices.

We also examined the level of congruence of gene partitions that are supported in each hypothesis across different analyses in Whelan2017_full matrices (Fig. Extended Data Fig. 3A). Consistent with the summary of $GLS analysis, We found gene population that supports in either hypothesis are largely congruent between site-homogenous and C models, whereas little overlap was observed in Poisson-CAT models (especially with different outgroup choice). In other words. our results strongly indicated that the distribution of phylogenetic signal of animal origin is largely sensitive to model selection (site-homogeneous vs. site-heterogeneous models). More importantly, the site-heterogeneous CAT model is also more sensitive to the outgroup choice than site-homogeneous and C60 models in Phillippe2009 and Whelan2017_full matrices.

Sensitive test of CAT categories {place holder for sensitive results}

The major difference between site-heterogeneous CAT and C60 models are the number of substitutional categories that are used in different data matrices. CAT models employ a Dirichlet process prior to infer number of substitutional categories or profiles with different nucleotide or amino acid frequencies to account for substitutional heterogeneity across a phylogenetic data matrix³ . Here, we reported the number of categories inferred from CAT model for each PhyloBayes analysis (Table XX). The number of categories is considerable larger than 60 categories used in C60 models, and in most cases even larger than the gene partitions within each data matrix (Table xx). Although the true number of substitutional categories in empirical data matrices is impossible to determine with current phylogenetic tools, we found that the amino frequencies of several categories inferred by C60 model are close to zero from IQtree results, potentially indicating that the C60 models might be already overfitting to most matrices. These results are consistent with the findings that CAT model often overestimate the true number of substitutional categories, especially for large data matrix from simulation analyses²⁴.

We next explore how the increase of number of categories may affect the phylogenetic inference related to the animal phylogeny with a panel of different substitutional categories ( in the representative data matrices (Fig. XX). Surprisingly, we found Poisson+nCAT60 showed exact same topology as Poisson+C60 or WAG+C60 with strong support of Ctenophora-sister, suggesting that the different phylogenetic analysis software (IQtree and PhyloBayes) and phylogenetic inference (ML and BI) has limited impact on phylogenetic inference to these matrices with similar evolutioanry models. Importantly, our results gave decreased support for Ctenophora-sister and increased support of Porifera-sister with the increase of number of substitutional categories (Fig. XX). More importantly, we are able to identify the transition of topology between the two alternative hypotheses in these three matrices. For example, the major transition occurred at 60 – 120 categories, and the support for Ctenophora-sister disappears after 120 categories in the Whelan2017_strict matrix. Interestingly, the results of our sensitive analyses further corroborated by the recent simulation study that the tendency of overestimate the true number of categories of Poisson-CAT or GTR-CAT model can potentially lead to a less accurate tree than the site-homogeneous models²⁴.

The GTR-CAT model often has a better model-fitting than Poisson-CAT model in previous analyses, and it is also suggested a better performance than Poisson-CAT model. Due to computational efficiency, we only ran GTR-CAT and GTR-nCAT60 models in the representative data matrices with the most closely related choanofalgellates as outgroups. Interestingly, we found both Whelan2017_full and Whelan2017_strict matrices shifted from Porifera-sister with Poisson-CAT model to Ctenophora-sister using GTR-CAT model. Moreover, we also found all results strongly supported the Ctenophora-sister with GTR-nCAT60 models. Thus, these results suggested that the Porifera-sister can only be recovered with several combined factors, including inclusive of a more simplistic evolutionary model, an unrealistic large number of substitution categories, more closely related outgroups and certain gene samplings, whereas Ctenophora-sister is supported in most other conditions.

Bias of data-recoding

There is growing interest in data-recoding and previous study⁹ hoped that data-recoding would reduce potential artifacts due to differences across species in amino acid frequencies. They report that posterior predictive (PP) analyses²⁷ indicate 6-state recoded analyses have better model adequacy than 20-state amino acid analyses, and “Porifera-sister was favored under all recoding strategies” in Whelan2015_D20 and Chang2015 data matrices. Here we focus on two aspects of Feuda et al. First, we found that many of their recoded analyses are actually unresolved (i.e., without strong support for either Porifera-sister or Ctenophora-sister), and that the analyses with the best posterior predictive scores do not provide strong support for Porifera-sister (Extended Data Fig. 3B). Second, we created four new random recoding schemes by shuffling the amino acids in the SR-6 scheme (see methods). When we applied each of these randomized codes to the Whelan matrix and analyzed them under the CAT-GTR+G model with phylobayes-mpi, we observed similar results as for the empirical SR-6 recoding. Like SR-6 recoding, random recoding increases support for Porifera-sister and improves the apparent adequacy of models to explain heterogeneity of states across taxa (PP taxon hetero mean and max, Extended Data Fig. 5). Thus, our analyses show the impact of recoding is largely due to discarding information, not accommodating variation in amino acid composition. These findings indicate that it is premature to accept Porifera-sister and reject Ctenophora-sister using data-recoding methods and recoding can be a problematic method for addressing compositional variation (see more detailed discussion in supplementary text).

Data selection and wrangling (Ben)

We retrieved matrices from each publication (Table 1), storing the raw data in this manuscript’s version control repository. We manually edited some minor formatting to make the batch processing of the matrices en masse, e.g. standardizing the formatting of charset blocks. All changes made are tracked with git.

Matrix comparison and annotation (Ben)

Phylogenetic analyses

MPKE AY IDGQ TRC WNLF SVH
EIFT WL QVPG NKM SCHA RYD
LCQS GK WPTI VRD YEFN MAH
IWQV TY KDLM ESH APCF GRN

Matrix mapping (Ben)

Taxa and partition correspondence across manuscripts was assessed by comparing all sequences for each taxon in each partition across all matrices with diamond blast. Based on inspection of sequence similarity, we excluded all comparisons with less than 99% identity and greater than 10^-25 e-value.

Taxa comparison across matrices (Ben)

The primary intent of comparing taxa across matrices was to validate our taxon name reconciliation across studies.

We first considered pairwise similarity between the same species from different matrices in different studies.

Models of molecular evolution

The exchangeability matrix \(E\) describes the rate at which one amino acid changes to another. Exchangeability matrices have been used in the studies under consideration here include: Poisson (or F81), WAG, LG, GTR. While the exchangeability matrix describes the relative rate of different changes between amino acids, the actual rate can be further scaled. There are a couple approaches that have been used in the studies considered here:

• F81³⁸ corresponds to equal rates between all states. The F81 matrix is also sometimes referred to as the Poisson matrix. It has no free parameters to estimate since all off-diagonal elements are set to 1.

• WAG³⁹ is an empirically derived exchangeability matrix based on a dataset of 182 globular protein families. It has no free parameters to estimate since all off-diagonal elements are set according to values estimated from this particular sample dataset.

• LG⁴⁰, like WAG, is an empirically derived exchangeability matrix. It is based on a much larger set of genes, and variation in rates across sites was taken into consideration when it was calculated. It has no free parameters to estimate since all off-diagonal elements are set according to values estimated from this particular sample dataset.

• GTR, the General Time Reversible exchangeability matrix, has free parameters for all off-diagonal elements that describe the exchangeability of different amino acids. It is constrained so that changes are reversible, i.e. the rates above the diagonal are the same as those below the diagonal. This leaves 190 parameters that must be estimated from the data long with the other model parameters and the phylogenetic tree topology. This estimation requires a considerable amount of data and computational power, but if successful has the advantage of being based on the dataset at hand rather than a different dataset (as for LG and WAG).

While the exchangeability matrix describes the relative rate of different changes between amino acids, the actual rate can be further scaled. There are a couple approaches that have been used in the studies considered here:

• Site homogeneous rates. The rates of evolution are assumed to be the same at all sites in the amino acid alignment.

• Gamma rate heterogeneity. Each site is assigned to a different rate class with its own rate value. This accommodates different rates of evolution across different sites. Gamma is used so commonly that sometimes it isn’t even specified, making it difficult at times to know if a study uses Gamma or not.

The vector of equilibrium frequencies \(\Pi\) describes the stationary frequency of amino acids. There are a few approaches that have been used across the studies considered here:

• Empirical site-homogeneous. The frequency of each amino acid is observed from the matrix under consideration and applied to homogeneously to all sites in the matrix.

• Estimated site-homogeneous. The frequency of each amino acid is inferred along with other model parameters, under the assumption that it is the same at all sites.

• CAT site heterogeneous³. Each site is assigned to a class with its own equilibrium frequencies. The number of classes, assignment of sites to classes, and equilibrium frequencies within the data are all estimated in a Bayesian framework.

• C10 to C60²². 10 to 60-profile mixture models as variants of the CAT model under the maximum-likelihood framework.

• Recoding method: Amino acids are recoded into six groups based on function to account for both compositional heterogeneity and substitution saturation. Several recoding strategies have been proposed, including Dayhoff 6-state recoding, S&R 6-state recoding, KGB 6-state recoding based on exchangeability matrix of models.

Details of published analyses

Data-recoding methods

Feuda et al.³⁶ tackled a difficult phylogenetic problem – placing the root of the animal phylogeny⁸. In the past decade, some analyses have placed Ctenophora (comb jellies) and others Porifera (sponges) as the sister group to all other animals (Fig. Extended Data Fig. 3B). Feuda et al.³⁶ present new phylogenetic analyses that they claim provide strong support for Porifera-sister. Their analyses consider datasets from studies by Chang et al.⁵³ and Whelan et al.⁵⁴, both of which found support for Ctenophora-sister. Feuda et al. were concerned that these Ctenophora-sister results were artefacts of lineage-specific differences in amino acid frequencies. In an attempt to reduce these differences, they recoded the full set of twenty amino acids into six groups of amino acids. These groups have more frequent evolutionary changes within them than between them, based on empirical observations in large protein datasets³⁵. The intent is to discard many lineage-specific changes, which are expected to fall within these groups. Rather than model compositional heterogeneity, as their title suggests, this approach discards heterogeneous information so that much simpler models with fewer states can be applied. They report that posterior predictive (PP) analyses²⁷ indicate 6-state recoded analyses have better model adequacy than 20-state amino acid analyses, and “Porifera-sister was favored under all recoding strategies”. Here we focus on two aspects of Feuda et al. First, we point out that many of their recoded analyses are actually unresolved (i.e., without strong support for either Porifera-sister or Ctenophora-sister), and that the analyses with the best posterior predictive scores do not provide strong support for Porifera-sister. Second, we present new analyses that show the impact of recoding is largely due to discarding information, not accommodating variation in amino acid composition. These findings indicate that it is premature to accept Porifera-sister and reject Ctenophora-sister. They also show that recoding can be a problematic method for addressing compositional variation.

Feuda et al. examine support for Ctenophora-sister and Porifera-sister under all combinations of two models of molecular evolution, four datasets, and four coding schemes. This provides 32 analyses that they report in their Table 3 and that we present graphically here as Extended Data Fig. 6. There is striking variation in support for Ctenophora-sister and Porifera-sister across these analyses (Extended Data Fig. 6). Feuda et al. accept the results of some analyses and reject others based on posterior predictive (PP) analyses of model adequacy, which assess how well a model explains variation in the data²⁷. They considered five different posterior predictive statistics that capture different types of variation in the data. From this they conclude that their “results strongly support sponges as the sister group of all other animals”.

This conclusion does not follow from their own presented results. Only a single analysis with posterior predictive scores provides what could be considered strong support > 95 posterior probability) for Porifera-sister. Of the 32 analyses, posterior predictive scores were calculated for 16 (those for the full Whelan and Chang matrices). Based on posterior predictive scores, Feuda et al. reject eight of these that were conducted under the GTR+G model (which all have strong support for Ctenophora-sister). This leaves eight CAT-GTR+G analyses (Extended Data Fig. 5). Two of these eight are analyses of the original 20-state amino acid data, both of which provide strong support for Ctenophora-sister. Of the six recoded analyses, five are unresolved. Only a single analysis for which posterior predictive scores are available provides strong support for Porifera-sister, the CAT-GTR+G analysis of the SR-6³⁵ recoded Whelan⁵⁴ matrix. Furthermore, this analysis does not have the best score according to any of the five posterior predictive statistics they considered (Extended Data Fig. 5). The only statistic that stands out for this one analysis is that it has the highest maxdiff (Extended Data Fig. 5), indicating that it did not converge as well as other analyses.

Though their study does not provide strong support for Porifera-sister, the sensitivity of their results to recoding provides an opportunity to better understand and evaluate the impact of recoding more generally. This is important given the growing interest in recoding³⁶. Feuda et al. hoped recoding would reduce potential artefacts due to differences across species in amino acid frequencies. They interpreted the fact that their analyses are sensitive to recoding as evidence that such an artefact exists and that they successfully addressed it by recoding. An alternative hypothesis is that recoding impacts phylogenetic analyses because it discards so much information. These two hypotheses can be tested by applying new recoding schemes that also reduce twenty states down to six, but are based on random grouping rather than empirical frequencies of amino acid exchange. Empirical and random recodings both discard the same amount of information, but only empirical recoding reduce the impact of amino-acid frequency as intended. Different results between empirical and random recoding would be consistent with the hypothesis that the empirical approach works as intended to accommodate compositional heterogeneity. Similar results would suggest that the impact of recoding is due to discarding information. Here we focus on their single analysis with a posterior predictive score that supports Porifera-sister, the CAT-GTR+G analysis of the SR-6 recoded Whelan data. we created four new random recoding schemes by shuffling the amino acids in the SR-6 scheme (see Supplemental Methods and analysis code at https://github.com/caseywdunn/feuda_2017). When we applied each of these randomized codes to the Whelan matrix and analyzed them under the CAT-GTR+G model with phylobayes-mpi³⁴, we observed similar results as for the empirical SR-6 recoding. Like SR-6 recoding, random recoding increases support for Porifera-sister and improves the apparent adequacy of models to explain heterogeneity of states across taxa (PP taxon hetero mean and max, Extended Data Fig. 3B).

These analyses suggest that the major impact of recoding on phylogenetic analyses is data reduction, not accommodation of compositional heterogeneity across species. This indicates that recoding may not be an effective tool for accommodating among-species differences in amino acid frequencies. Compositional heterogeneity would be better addressed with models of molecular evolution that explicitly describe such differences⁵⁵, if progress can be made on the considerable computational challenges of such complex models.

Partition comparison across matrices (XXBen, do we need this?)

The count for a partition pair can be much larger than the number of genes in the matrix, which suggests that the count is the number of hsps rather than the number of sequences with hits.

There are 17 matrices. A gene that is perfectly sampled would form a cluster with this size. To our surprise, very few clusters, though, are this size. This suggests that intersection of genes between matrices is extreme low. This finding suggested genes used for exploring the position of the root of the animal phylogeny